
clear all;
clc;
disp('In-sample optimization to find w = arc min f(X,Z), where f is a vectors. We use one of multi objective programming approach in Matlab fminimax');
load dataAnalysisQ_G1000k.mat

INSAMPLE_YR = 2007;
INSAMPLE_QR = 4;
CRITERIA = 11;

inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = length(SVD)-inSampleLength;

xx =zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
avgRollingTM =zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
bin=zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
InsampleSVD = zeros(inSampleLength,1);
InsampleMobilityNorm = zeros(inSampleLength,1);
PHat=zeros(10,inSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
N = zeros(inSampleLength+outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
naive1 = zeros(length(SVD), RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );
naive2 = zeros(length(SVD), RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );

filename='optimization_08q4_jun13.xls';


%To calculate the rolling avg TM and rolling credit score/bins-------------

for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:,:) = durationCount(i,j,:,:);       
        xx = xx + x;
        t = (i-1)*4 + j;
        N(t,:,:)= x;
        
        for ii = 1 : RATE_LIST_LENGTH-1
            for jj = 1: RATE_LIST_LENGTH
                avgRollingTM (t, ii, :) = xx(ii,:)./ sum(xx(ii,:),2);
                cdf = sum( avgRollingTM(t, ii, jj:RATE_LIST_LENGTH) );
                bin(t,ii,jj)=norminv(cdf,0,1);
            end
        end
        
    end
end

disp('finish Rolling Avg Transition Matrix and Bin');


% To calculate credit cycle index Z, real and forecast----------------------
% We here use empirical Z to calculate w.

% Take out the in-sample data

% speculative 3.5-4.5; investment 1.5-3
defaultFreqInv = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);
defaultFreqSpe = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);


for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:, :) = durationCount(i,j,:,:);
        defaultFreqInv(4*(i-1)+j,:) = sum(x(1:4,RATE_LIST_LENGTH)) / sum(sum(x(1:4,:)));
        defaultFreqSpe(4*(i-1)+j,:) = sum(x(5:7,RATE_LIST_LENGTH)) / sum(sum(x(5:7,:)));
    end
end


PD = defaultFreq(1:inSampleLength);
PDInv = defaultFreqInv(1:inSampleLength);
PDInv(2,:) = 0.00001; % the empirical value is 0, so we artifically change it to extremely small number
PDSpe = defaultFreqSpe(1:inSampleLength);
inversePD = norminv(PD);
inversePDInv = norminv(PDInv);
inversePDSpe = norminv(PDSpe);

InsampleAvgTM (:,:) = avgRollingTM (inSampleLength, :, :);
InsampleAvgTMmobilityNorm = mean(svd(InsampleAvgTM));

for i = 1: inSampleLength
    b(:,:) = stdProbQ (i, :, :);
    InsampleMobilityNorm (i,1) = mean(svd(b));
    InsampleSVD(i,1) = InsampleMobilityNorm (i,1) - InsampleAvgTMmobilityNorm; %SVD
end

%impirical Z value

Z_PD = - zscore(inversePD); 
Z_SVD = - zscore(InsampleSVD);



% To minimize the distance of each time point to get estimate of w------------
% PHat_BFS = zeros (inSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);

R = RATE_LIST_LENGTH-1;
C = RATE_LIST_LENGTH;
P = stdProbQ;
X(:,:) = bin(inSampleLength,:,:);
% wantedZ = Z_PD; 
wantedZ = Z_SVD;
t = inSampleLength;
W = zeros(CRITERIA,1);

lb = zeros(1,inSampleLength);
ub = ones(1,inSampleLength);
A = []; b =[];
Aeq =[];beq=[];
x0 = 0.2;

% L1
obj1=@(w)difL1_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(1,1),fval] = fminimax(obj1,x0,A,b,Aeq,beq,lb,ub);


% L2
obj2=@(w)difL2_MultiObj(w,t,N,X,wantedZ,R,C); 
[W(2,1) ,fval] = fminimax(obj2,x0,A,b,Aeq,beq,lb,ub);


% WAD
obj3=@(w)difWAD_MultiObj(w,t,P,X,wantedZ,R,C); % weight is est. P
[W(3,1),fval] = fminimax(obj3,x0,A,b,Aeq,beq,lb,ub);


% NAD
obj4=@(w)difNAD_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(4,1),fval] = fminimax(obj4,x0,A,b,Aeq,beq,lb,ub);


% SVD
obj5=@(w)difSVD_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(5,1),fval] = fminimax(obj5,x0,A,b,Aeq,beq,lb,ub);


% NSD distance
obj6=@(w)difNSD_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(6,1),fval] = fminimax(obj6,x0,A,b,Aeq,beq,lb,ub);


% BFS distance
obj7=@(w)difBFS_MultiObj(w,t,N,P,X,wantedZ,R,C); 
[W(7,1) ,fval] = fminimax(obj7,x0,A,b,Aeq,beq,lb,ub);


 % D1^2
obj8=@(w)difD1sqr_MultiObj(w,t,P,X,wantedZ,R,C); % use 2 Z index
[W(8,1),fval] = fminimax(obj8,x0,A,b,Aeq,beq,lb,ub);


% D2^2
obj9=@(w)difD2sqr_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(9,1),fval] = fminimax(obj9,x0,A,b,Aeq,beq,lb,ub);


% D3^2
obj10=@(w)difD3sqr_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(10,1),fval] = fminimax(obj10,x0,A,b,Aeq,beq,lb,ub);

% D4^2
obj11=@(w)difD4sqr_MultiObj(w,t,P,X,wantedZ,R,C); 
[W(11,1),fval] = fminimax(obj11,x0,A,b,Aeq,beq,lb,ub);



disp('finish w estimation');

% To calculate in-sample fitted matrix and error(Mean Absolute Error)-------------------

for i= 1: CRITERIA
    
  PHat(i,:,:,:) = calculatePHat_fixedW(W(i,1),t,X,wantedZ,R,C); % PHat = f(X,What,Z)
    
end


for i = 1: NUMBER_OF_YEAR*4
if i == 1
    naive1(i,:,:) = avgRollingTM (i,:,:); % take average as naive1 benchmark
else
    naive1(i,:,:) = avgRollingTM (i-1,:,:);
end
end

for i = 1: NUMBER_OF_YEAR*4
if i == 1
    naive2(i,:,:) = stdProbQ (i,:,:);% take previous as naive2 benchmark
else
    naive2(i,:,:) = stdProbQ (i-1,:,:);
end
end

e = zeros(CRITERIA,inSampleLength,3);
Error = zeros(CRITERIA,3);
avgPDfcst = zeros(CRITERIA,t,4); % avg default prob for model, benchmarks and empirical

for i = 1:CRITERIA
    
    for j = 1:t
    num(:,:)= N(j,:,:);   
    x1(:,:) = PHat(i,j,:,:);
    x2(:,:)= naive1(j,:,:);
    x3(:,:) = naive2(j,:,:);
    y(:,:) = stdProbQ(j,:,:);
    
    e(i,j,:) = calculateMAE(i,x1,x2,x3,y,R,C);
    
     avgPDfcst (i,j,1) = sum(x1(:,C).* sum(num,2))./sum(sum(num)) ;
     avgPDfcst (i,j,2) = sum(x2(:,C).* sum(num,2))./sum(sum(num)) ;
     avgPDfcst (i,j,3) = sum(x3(:,C).* sum(num,2))./sum(sum(num)) ;
     avgPDfcst (i,j,4) = sum(y(:,C).* sum(num,2))./sum(sum(num)) ;
      
    end
    
    Error (i, 1) = mean (abs (e(i,2:inSampleLength,1)));
    Error (i, 2) = mean (abs (e(i,2:inSampleLength,2)));
    Error (i, 3) = mean (abs (e(i,2:inSampleLength,3)));
    
end
disp('finish error calc');

% visual check est. TM-----------------------------------------------------

% first, we look at cell prob. 

cellProb5 = cellProbCompare(PHat,naive1, naive2, stdProbQ, R, C, inSampleLength,5, 0);
cellProb9 = cellProbCompare(PHat,naive1, naive2, stdProbQ, R, C, inSampleLength,9, 0);
cellProb10 = cellProbCompare(PHat,naive1, naive2, stdProbQ, R, C, inSampleLength,10, 0);


linestyles = cellstr(char('-',':','-.','--','-',':','-.','--','-',':','-',':',...
'-.','--','-',':','-.','--','-',':','-.'));
Markers=['o','x','+','*','s','d','v','^','<','>','p','h','.',...
'+','*','o','x','^','<','h','.','>','p','s','d','v',...
'o','x','+','*','s','d','v','^','<','>','p','h','.'];
MarkerEdgeColors=jet(5);  % n is the number of different items you have

% Avg fcst PD graph 
figure
f(:,:) = avgPDfcst(5,:,1);
plot (f,[linestyles{1} Markers(1)]);
legend('model w/ SVD','Location','NorthWest');
set(gca,'XTick',1:4:inSampleLength);
labels = quaterlabels(1998, inSampleLength);
labels = labels(1:4:inSampleLength);
set(gca,'XTickLabel',labels);
xlabel('Quarter');
title('Average Rating Default Probability Estimation: 1998Q1 - 2007Q4');
hold on
for i = 1:4
plot(avgPDfcst(7+i,:,1),[linestyles{1} Markers(i+1)]);
end
hold on
for i = 2:4
plot (avgPDfcst(5,:,i),[linestyles{i}],'Color',MarkerEdgeColors(i,:));
end

% SVD criteria Default Prob
figure
for i = 1:R
temp(:,:) = cellProb5(i,C,:,:);
subplot(4,2,i); plot (temp);
% title('In-sample SVD fitted PD Comparision, Rating',... 
%   'FontWeight','bold')
set(gca,'XTick',1:4:inSampleLength);
labels = quaterlabels(1998, inSampleLength);
labels = labels(1:4:inSampleLength);
set(gca,'XTickLabel',labels);
% xlabel('Quarter');
end
title('In-sample SVD fitted PD Comparision, Rating',... 
  'FontWeight','bold')
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

% SVD criteria diagonal Prob
figure
for i = 1:R
% figure;
temp(:,:) = cellProb5(i,i,:,:);
subplot(4,2,i); 
plot (temp);
% title('In-sample SVD fitted Diagonal Probability Comparision',... 
%   'FontWeight','bold')
set(gca,'XTick',1:4:inSampleLength);
labels = quaterlabels(1998, inSampleLength);
labels = labels(1:4:inSampleLength);
set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
end
title('In-sample SVD fitted Diagonal Probability Comparision',... 
  'FontWeight','bold')
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');


% D2sqr criteria Default Prob
figure
for i = 1:R
temp(:,:) = cellProb9(i,C,:,:);
subplot(4,2,i); plot (temp);
title('In-sample D2 fitted Default Probability Comparision',... 
  'FontWeight','bold')
set(gca,'XTick',1:4:inSampleLength);
labels = quaterlabels(1998, inSampleLength);
labels = labels(1:4:inSampleLength);
set(gca,'XTickLabel',labels);
% xlabel('Quarter');
end
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

% D2^2 criteria diagonal Prob
figure
for i = 1:R
% figure;
temp(:,:) = cellProb9(i,i,:,:);
subplot(4,2,i); 
plot (temp);
% title('In-sample D2 fitted Diagonal Probability Comparision',... 
%   'FontWeight','bold')
set(gca,'XTick',1:4:inSampleLength);
labels = quaterlabels(1998, inSampleLength);
labels = labels(1:4:inSampleLength);
set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
end
title('In-sample D2 fitted Diagonal Probability Comparision',... 
  'FontWeight','bold')
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');


% % D2sqr criteria column rating 4
% figure
% for i = 1:R
% temp(:,:) = cellProb9(i,4,:,:);
% subplot(4,2,i); plot (temp);
% title('In-sample D3 fitted Rating 4 Comparision for rating',... 
%   'FontWeight','bold')
% set(gca,'XTick',1:4:inSampleLength);
% labels = quaterlabels(1998, inSampleLength);
% labels = labels(1:4:inSampleLength);
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% end
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');














% % criteria D3
% g1 (:,1)= PHat(8,:,1,1); g1(:,2) = naive1(:,1,1); g1(:,3) = naive2(:,1,1); g1(:,4)= stdProbQ (1:length(PHat),1,1);% rank 1.5 - 1.5 
% g2 (:,1)= PHat(8,:,3,3); g2(:,2) = naive1(:,3,3); g2(:,3) = naive2(:,3,3); g2(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5
% g3 (:,1)= PHat(8,:,6,6); g3(:,2) = naive1(:,6,6); g3(:,3) = naive2(:,6,6); g3(:,4)= stdProbQ (1:length(PHat),6,6);% rank 3 - 3
% g4 (:,1)= PHat(8,:,5,8); g4(:,2) = naive1(:,5,8); g4(:,3) = naive2(:,5,8); g4(:,4)= stdProbQ (1:length(PHat),5,8);% rank 3.5 -D
% g5 (:,1)= PHat(8,:,6,8); g5(:,2) = naive1(:,6,8); g5(:,3) = naive2(:,6,8); g5(:,4)= stdProbQ (1:length(PHat),6,8);% rank 4 -D
% g6 (:,1)= PHat(8,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
% 
% % criteria SVD
% 
% g6_svd (:,1)= PHat(5,:,7,8); g6_svd(:,2) = naive1(:,7,8); g6_svd(:,3) = naive2(:,7,8); g6_svd(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
% g2_svd (:,1)= PHat(5,:,3,3); g2_svd(:,2) = naive1(:,3,3); g2_svd(:,3) = naive2(:,3,3); g2_svd(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5
% 
% s =  zeros(t,4);
% 
% for i = 1:t
%     x1(:,:) = PHat(5,i,:,:); x2(:,:) = naive1(i,:,:); x3(:,:) = naive2(i,:,:); y(:,:) = stdProbQ (i,:,:);
%     s(i,1) = mean (svd(x1));
%     s(i,2) = mean (svd(x2));
%     s(i,3) = mean (svd(x3));
%     s(i,4) = mean (svd(y));
% end
% 
% figure
% plot (g6); %([g1,gg1,ggg1]);%,'-*'); plot (gg1,'-rs'); 
% set(gca,'XTick',1:4:length(g1));
% labels = quaterlabels(1998, length(g1));
% labels = labels(1:4:length(g1));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

% then we look at metrics

      
% % D3
% obj8=@(w)difD3_MultiObj(w,t,P,X,wantedZ,R,C); % use 2 Z index
% [W(8,1),fval] = fminimax(obj8,x0,A,b,Aeq,beq,lb,ub);
% 
% 
% % D1
% obj9=@(w)difD1_MultiObj(w,t,P,X,wantedZ,R,C); 
% [W(9,1),fval] = fminimax(obj9,x0,A,b,Aeq,beq,lb,ub);
% 
% 
% % D1sqr
% obj10=@(w)difD1sqr_MultiObj(w,t,P,X,wantedZ,R,C); 
% [W(10,1),fval] = fminimax(obj10,x0,A,b,Aeq,beq,lb,ub);

 
